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ISU — Multigrid for computing propagators 
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The Iteratively Smoothing Unigrid algorithm (ISU), a new multigrid method for computing propagators in 
Lattice Gauge Theory, is explained. The main idea is to compute good (i.e. smooth) interpolation operators in 
an iterative way. This method shows no critical slowing down for the 2-dimensional Laplace equation in an SU(2) 
gauge field. First results for the Dirac-operator are also shown. 



1. The problem 

The greatest obstacle for doing realistic (un- 
quenched) simulations of lattice QCD is the 
large computer time required for the inversion of 
the Dirac-operator due to critical slowing down. 
Multigrid algorithms, which have been successfull 
for the solution of differential equations describ- 
ing ordered systems, have been studied for some 
time, but only with limited success. 

In jj| we proposed the ideas of a new multi- 
grid method to overcome these problems. This 
new algorithm, called Iteratively smoothing uni- 
grid or ISU has been described in detail in M. 
Here we want to review it shortly and on general 
grounds and want to investigate its performance 
for the case of the 2-dimensional SU(2) bosonic 
and fcrmionic propagator equation. 

2. The Unigrid 

In this section, the basic principles of the multi- 
grid method are presented from the unigrid point 
of view. 

The two key observations are: 

1. Standard relaxation algorithms (like Gaufi- 
Seidel-relaxation) applied on a lattice with 
lattice constant a are not efficient in reduc- 
ing, but in smoothing the error on this scale. 
(The error is defined as the difference be- 
tween the approximate and the exact solu- 
tion.) 



2. A function which is smooth on a length 
scale a can be obtained by interpolation 
from a grid with lattice constant ex a. 

So the idea is to introduce, in addition to the 
fundamental lattice A on which the problem is 
defined, auxiliary layers A J with lattice constants 
dj = L J b ao, where Lf,, the blocking factor, is 2 or 
3. The last of these lattices, A N , consists of only 
one point. Points x S A J can be identified with 
the corresponding points in A . We define blocks 
[x] in A with sidelength 2a j — 1 and the point x 
in the center. 

Let the fundamental equation be D£ = f . 

To an approximate solution £ corresponds an 
error e = £ — £. We can cast the equation into the 
form D e = r, where r = f — D £ is the residual. 

After relaxing on A , e is smooth on scale ao. It 
can therefore be obtained by interpolating it from 
a function e 1 on A 1 with the help of the interpo- 
lation operator A^° 1 ' . ^4'° J ' maps functions on A- 7 
to functions on A , its adjoint A^°^ blocks func- 
tions from A to A- 7 . It has the property A zx = 
for z £ [x]. So we get e 
the coarse grid equation 

_ 4 [oi]* D ^[oi] e = ^[oi]* ] 



^[oi] gi^ resulting in 
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On this equation we relax again, thereby 
smoothing the error on scale a±. We correct the 
approximation on the fundamental lattice £ <-^ 
£ — A^° ^ e 1 and then we go to the next-coarser 
layer. (This is the difference to a true multigrid 
method: There one goes from A 1 directly to the 
next layer, without correcting first on the funda- 



mental lattice. This is advantageous, because the 
work involved is less, but that the ISU-algorithm 
cannot work as a true multigrid, see |2j.) 

3. The need for smooth operators 

The method described above can only work 
well if we know what is meant by the smooth- 
ness of the error after relaxation. It can be easily 
seen that in the context of lattice gauge theo- 
ries, smoothness is not a priori defined, because 
no unique way to compare function values at dif- 
ferent sites can be defined. Therefore we have 
to look for a new, appropriate definition. The 
following two definitions are in most cases equiv- 
alent. 

Def. 1: A mode £ is smooth on scale A iff || L£ || <C 
£, L in units A = 1 

Def. 2: A mode is smooth on scale A iff it is 
not efficiently reduced by relaxation on this 
scale. 

That relaxation produces a smooth error not 
only according to the second, but also to the first 
definition, can be seen from the fact that usu- 
ally the low-lying eigenmodes of D are the bad- 
converging ones. (See H for a discussion of this 
and some caveats.) 

To interpolate an error smooth in this sense, 
the interpolation operators A^°^ have to be 
smooth themselves. So we have to find out how 
to calculate smooth interpolation operators. As 
the smoothest mode (according to def. 1) is the 
lowest eigenmode of D, it seems that we have to 
solve an equation as difficult as the one we started 
with to find the .A^. 

4. The iteratively smoothing unigrid 

But this is not true because of the following 
observation: 

It is easier to calculate the shape of a mode 
which converges badly in a given iteration scheme 
than to reduce it directly with this scheme. 

This can be seen from the following exam- 
ple: Imagine that D only has one bad-converging 
mode ijj. Then we can compute this easily by 



trying to solve D£ = 0, because in this case the 
approximate solution equals the error, which con- 
verges fast to a (small) multiple of ip. 

The second problem solving strategy is to solve 
our problem by reducing it to many similar prob- 
lems which can be solved step by step. 

The ISU algorithm for calculating good inter- 
polation operators A^- -*' works as follows: 

To calculate A J on a support [a;]: 

Solve D|r x ]„4L; = SqA z , j I via inverse iteration, 
using all interpolation operators on layers A fe with 
k < j. (The righthandside could also be taken to 
be zero as in the above example, according to def. 
2 of smoothness.) I^j here means restriction to 
the block [x], using e.g. Dirichlet boundary con- 
ditions. 

The crucial point here is that the already cal- 
culated operators on the finer layers are used for 
the iteration, otherwise there would be many bad- 
converging modes and the calculation of A*- -*' 
would be slow. 

5. Performance of ISU 

We studied the performance of this algorithm 
in a two-dimensional SU(2) lattice gauge field, us- 
ing the equation 



(D - £ + Sm 2 )£ = f 



(2) 



where D was chosen to be the negative Laplace 
(bosonic case) or the negative squared staggered 
Dirac operator. £o is the lowest eigenvalue of the 
operator which was subtracted to be able to tune 
criticality by changing the parameter Sm 2 . This 
subtraction is necessary for the Laplace operator 
(because its lowest eigenvalue can be large) and 
it eases the analysis of the critical behaviour also 
for the Dirac case. 

We measured the asymptotic convergence rate 
t, i.e. the number of iterations needed to reduce 
the error by a factor e, where by an iteration we 
mean visiting each layer of the unigrid twice in a 
so-called V(l,l)-cycle (see (§). 

In the bosonic case we studied the convergence 
rate on grids of sizes 32 2 -128 2 at various values 
of the inverse coupling /3. It was found, that t 
was approximately 1, independent of the lattice 
size, /3, and Sm 2 (for small enough Sm 2 ) which 



Table 1 

Convergence rate r of the ISU-algorithm applied 
to the squared Dirac equation as a function of the 
grid length for physically scaled gauge fields. 



Grid size 18 2 



54 2 



162 2 



# runs 


50 


50 


9 





1 


9 


81 


Sm 2 


0.002 


0.002/9 


0.002/81 


T 


25.6 ±0.9 


6.67 ±0.14 


7.5 ±0.8 



means that the algorithm eliminates critical slow- 
ing down completely. Moreover it was found that 
the number of inverse iterations to calculate the 
smooth interpolation operators was six, again in- 
dependent of the problem parameters, so critical 
slowing down does not enter through the back- 
door. A detailed discussion of this result can be 
found in g. 

In applying the algorithm to fermions we used 
a block-factor of 3 instead of the more usual 2 to 
preserve the symmetry of the staggered grid also 
on the coarser layers H . The grid length therefore 
should be chosen as L — 2 ■ 3 N , where N is the 
number of layers. To improve convergence, we 
introduced an additional layer A add , consisting of 
only four points, one for each pseudoflavour. This 
allows for more interpolation operators that cover 
a large part of the grid. (In the case of bosons, one 
could also use a block factor of three. There, the 
additional layer is not needed, because we found 
that even for this block factor r was quite small 
and there was no critical slowing down.) 

Table [l] shows preliminary results for the con- 
vergence rate for physically scaled gauge fields, 
i.e. for fields where (3 ex L 2 and 5m 2 ex L~ 2 . 
The calculations were done on grids with grid size 
18 2 -162 2 with P = 1.0 and Sm 2 = 2 • 10~ 3 on 
the smallest lattice. Obviously there is no critical 
slowing down in this sense, which should be ex- 
pected as the squared Dirac operator approaches 
the Laplacian in the continuum limit. 

However, for small j3 the convergence rates are 
bad (large absolute values) even on small lattices. 
The situation is much worse than for bosons, 



where r was practically independent of all pa- 
rameters. 

6. Conclusions 

The great success of the algorithm in case of 
the Laplace equation and the results obtained for 
the Dirac case pose three questions: 

1. What is the crucial difference between the 
Laplacian and the Dirac operator, causing 
the one to converge much better than the 
other? 

2. Can the algorithm be improved, perhaps by 
including insights gained from the answer to 
the first question? 

3. How does the algorithm (or its improved 
version) behave in four dimensions? 

At the moment we are studying mainly the first 
question. It seems that at least part of the answer 
lies in the shape of the low-lying eigenmodes of 
the Dirac-operator. 
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